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Abstract A new class of third order Runge-Kutta methods for stochastic differential 
equations with additive noise is introduced. In contrast to Platen's method, which to 
the knowledge of the author has been up to now the only known third order Runge- 
Kutta scheme for weak approximation, the new class of methods affords less random 
variable evaluations and is also applicable to SDEs with multidimensional noise. Or- 
der conditions up to order three are calculated and coefficients of a four stage third 
order method are given. This method has deterministic order four and minimized 
error constants, and needs in addition less function evaluations than the method of 
Platen. Applied to some examples, the new method is compared numerically with 
Platen's method and some well known second order methods and yields very promis- 
ing results. 
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1 Introduction 

In many applications, e. g., in epidemiology and financial mathematics, taking stochas- 
tic effects into account when modelling dynamical systems often leads to stochastic 
differential equations (SDEs). An important subclass of these are SDEs with additive 
noise in the form 

rt m 

X(t)=x + go(s,X(s))ds + Tg { (Fi(t)-Wi(tQ)). (1.1) 
A) ,=i 
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Here, W(t) is an m-dimensional Wiener process defined on a probability space 
(£2, A, the Borel-measurable drift go : W* — > W 1 is assumed to be sufficiently 
differentiable and to satisfy a Lipschitz and a linear growth condition, and gi € M. d , i = 
l,...,m. Then the Existence and Uniqueness Theorem [ 10 1 applies. Examples of such 
systems arising in experimental psychology, turbulent diffusion, radio-astronomy and 
blood clotting dynamics can be found in ifTTl . 

In recent years, the development of numerical methods for the approximation 
of SDEs has become a field of increasing interest, see e. g. Ill 111151 and references 
therein. Whereas strong approximation methods are designed to obtain good pathwise 
solutions, see e. g. J3), weak approximation focuses on the expectation of functionals 
of the solution: 

Let Cp(R d ,R) denote the space of all g E C l (R d , R) fulfilling a polynomial growth 

condition [11). Further, let I 1 ' = {?o,?i, ,?at} with to < t\ < . . . < f# = T be a 

discretization of the time interval / = [to,T] with step sizes h„ =t n +\ — t„ for n = 
0,l,...,iV-l. 

Definition 1 (weak convergence) A time discrete approximation Y h — (Y h (t)) teI h 

converges weakly with order p toX ash-^0 at time t £ I h if for each f G Cp + (M d ,R) 
there exist a constant Cf and a finite So > such that 



holds for each h €]0, 8q[. 

Many approximation schemes for SDEs fall into the class of stochastic Runge- 
Kutta (SRK) methods. Second order SRK methods for the weak approximation of 
SDEs were proposed by Kloeden and Platen IfTTl . Komori T\3l . Mackevicius and 
Navikas QU, Tocino and Vigo-Aguiar ED, RoBler II18II191 . and Debrabant and RoGler 
CIEIIH. An explicit third order weak SRK method for autonomous SDEs with addi- 
tive scalar noise as well as its generalization to general scalar noise have been given 
in Kloeden and Platen IfTTl . However, the authors state there that "it remains an open 
and challenging task to derive simpler derivative free order 3.0 weak schemes, at least 
for important classes of stochastic differential equations." The present article solves 
this problem in the case of additive noise and overcomes also the restriction to scalar 
additive noise. 

To do so, we consider the following class of s-stage SRK methods, 



which defines a li-dimensional approximation process Y h with Y h (t„) = Y„. Here, 
Jk, k — 1 , . . . , 2m, are independent random variables which do not depend on h„ 
and whose moments all exist. Further, a = (ai , . . . , Ot s ) J , A — (fli/)f /=i c — 
(ci,...,Cj) T , b\ = (fei,i,...,Ai,i) T , and b 2 = (Z>2,i, ■ ■ ■ M.s) J are the coefficients of 



\E(f(Y h (t)))-E(f(X(t)))\<C f hP 




(1.2a) 




(1.2b) 
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the SRK method. In the following we choose c = All with H = (1, . . . , 1) T G W. Con- 
sequently, from now on we can assume for the analysis of this methods that SDE 
( 11.11 is given in autonomous form, i. e., go(t,X) = go(X). The analysis relies on the 
theory of stochastic B-series, which is shortly reviewed in Section [2] and applied in 
Section|3]to derive order conditions for method (11.2b up to order three. Then, in Sec- 
tion |4] a concrete explicit third order method is constructed by minimizing the error 
coefficients. Finally, in Section|5]we give some numerical examples. 

2 Stochastic B-series 

Order conditions for method ( 11.21 ) can be calculated using the colored rooted tree 
theories derived for the weak approximation of Ito respectively Stratonovich SDEs by 
SRK methods, compare ni6U17lfT2l . Here, we will follow the more general approach 
developed in [5 1, which is based on the work in [1.2, 17] and applicable both for Ito- 
and Stratonovich SDEs as well as strong and weak approximation. For more details 
and proofs, see 0. 

First, we introduce the set of colored, rooted trees related to the SDE ( 11.11 ). as 
well as the elementary differentials associated with each of these trees. We adapt 
these definitions to the special case of additive noise by neglecting all terms which 
are related to derivatives of gi, I = 1 , . . . ,m. 

Definition 2 (trees) The set ofm + l-colored, rooted trees 
T add = {<d}UT U {•!,...,.,„} 

related to additive noise is recursively defined as follows: 

(a) The graph »o = [0]o with only one vertex of color belongs to Tq. 

Let T = [ti, T2, . . . , Tk-]o be the tree formed by joining the subtrees T\, T2, . . . ,T K each 
by a single branch to a common root of color 0. 

(b) //Ti , T 2 , • • • , T K € T add , then T = fa , T 2 , • • • , T K ]o € T . 

Thus, Tq is the set of trees with a 0-colored root. »o will be called deterministic node, 
•/ for I > stochastic node of color I. 

Definition 3 (elementary differentials) For a tree % e T add the elementary differen- 
tial is a mapping F{%) : W 1 — > W ! defined recursively by 

(a) F(9)(xq)=xq, 

(b) F(»o)(xo) =go(xo), F(9i)(x )=giforl = l,...,m, 

(c) IfT = [zi,x 2 ,---,*K]oeT ,then 

F(T)( Xo )^g^\x Q )(F(T l )(x ),F(T 2 )(x Q ),...,F(T K )(x )). 

To simplify the presentation, we neglect in the following the index n of h n and write 
only h. Further, we denote by S the set of families of Borel measurable mappings 

S := {{(p(h)} h >o ■ <p{h) ■ ^ Kis ^-^-measurable Vh > 0}. 

Both the solution of ( II. lb and its approximation by method ( 11.2b can formally be 
written in terms of B-series. 
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Definition 4 (B-series) Given a mapping (j) : T — > S satisfying 
0(0) = 1 and <j)(x)(0) =0, Vt G T add \{®}. 
A (stochastic) B-series is then a formal series of the form 

B(>,*O;/0= L OC{x)^(x)(h)-F(x)(x ), 

zeT add 

where a : T — > Q is given by 

1 K 

o(fl) = l, (*(•/) = 1, o(t= [Ti,...,T r ]i) = — — rn a ( T >)' 

where ri,rz, ■ ■ ■ ,r q count equal trees among Ti, T2, . . . ,X K . 
For multidimensional : T" dd —> E s , s G N, we define 

B{<t>,x ;h) = [B(^ u x ;h),...,B(^ s ,x ;h)] T . 

If Z(n) can be written as a B-series, then f(Z(h)) can be written as a similar series, 
where the sum is taken over trees with a root of color / and subtrees in T add : 

Lemma 1 IfZ(h) = B(<j>,x ;h) is some B-series and f G C°°(W',R d ), then f{Z(h)) 
can be written as a formal series of the form 

f(Z(h)) = £ p(u)-w(u)(h)-G(u)(x ), (2.1) 

ueuf d 

where 

(a) Uf dd is a set of trees derived from T add as follows: [0]/ G U add , and ifx\ , Xz, . . . , X K G 
T add , then [ti , X 2 , . . . , X K ] f G U a f dd , 

(b) G([%)(io)=/(io)flndG(^[T 1 ,...,T r ] / )(io)-/ (,c) (io)(f(Ti)(^),-,f(^)(^)), 

(c) j3([0]/) = 1 and P(u = [x 1 ,...,X K ] f ) = n! ,- 2 l... r<?! nj = i«(^)» w/iere ri,r 2 ,..., 
r q count equal trees among X\,X%,...,X K , 

(d) V^([0]/) = 1 and y^u = [x h . . . , x K ] f )(h) = njLi (*)• 

Remark 2.1 To simplify the presentation, we assume throughout this article that all 
derivatives of / and go exist. Otherwise, one had to consider truncated B-series with 
a remainder term. 

Theorem 2 The solution X(tQ + h) of ( 11. Il l can be written as a B-series B(q>,xo;h) 
with 

<p(0)=l, <p(. Q )(h)=h, <p(.,)(h)=W,(h), l = l,...,m, 

rh K 

(p(x = [x 1 ,...,x lc ] )(h)= / Y\(p(Xj)(s) ds. 

JO j=l 

The following definition of the order of the tree, p(x), is motivated by the fact that 
EW,{h) 2 =hforl> 1. 
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Definition 5 (order) The order of a tree % £ T add is defined by 
p(0) = O, P{*i) = \, l = l,...,m, 

and 

' 1 for 1 = 0, 



P(t = [ti,...,t k ] / ) = £p(t ! -) 



i (J /or/>0. 

r/je order of a tree u £ Uf dd is given by p(u = [ti, . . . , T K ] / ) = £ p (t,). 

In the following we define the product of vectors by componentwise multiplica- 
tion. 

Theorem 3 The numerical approximation Y\ as well as the stage values can be writ- 
ten in terms ofB-series 

H = B(® H ,x ;h), Yi=B{®,x ;h) 

with 

4>tf(0) = l, <P H {» l )(h) = Vh{b 1 J l + b 2 J m+ i), l=l,...,m, (2.2a) 

%(t=[ tlr .,i r ]o)(A)=MnJ =1 %(T/)(li) (2.2b) 

and 

0(0) = 1, <P(*i)(h) = VhJ,, l = l,,,,,m, (2.3a) 

*(t = [x x ,.. .,T K ] )(h) = ha T UU ®n(Tj)(h). (2.3b) 



3 Derivation of order conditions 

With all the B-series in place, we can now present the order conditions for the weak 
convergence. 

Let lef(h;t,x) be the weak local error of the method starting at the point (t,x) 
with respect to the functional / and step size h, i. e. 

le f {h;t,x) =E(f(Y h (t + h))-f(X(t + h))\Y h {t) =X(t)=x). 

From Theorems [2] and |3] and Lemma[T]we obtain 

le f (h;t,x)= £ P(u)-B[Y 9 (u)(h)-Y 9 {u){h)] -G(u)(x) (3.1) 

ueuf d 

with 

Yr,([8]/) = 1, W9(u = [ri,...,r K ]f)(h) = l\( P (Tj)(h) (3.2) 

and 

^([%) = 1, W<t ,(u=lT u ...,T K } f )(h) = fl<P(T J )(h). (3.3) 
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Thus, we have weak consistency of order p (and thus, due to the Milstein theorem 
03), also weak convergence) if and only if 

E\i/4,(u)(h)=E\i/Ju)(h) + @(h p+l ) \fu£U l ) dd with p{u) <p + \. (3.4) 

2 

Note that (13.4b slightly weakens conditions given in ifTTIl . 

By Theorems |2] and [3] d3.21 > and (13.31 > we can now evaluate the order conditions 
d3.41 > and obtain the following theorem. 

Theorem 4 For a p-th order method choose the independent random variables 
of the SRK method (11.21 ) such that their moments coincide with those of N {0,1) 
up to the (2p + l)-th moment for k=l,. . .,m and up to the [2p — l)-th moment for 
k=m+l,. . . ,2m. If in addition the coefficients of the SRK method (11.2b fulfill 

1. a J l = 1, 

then the method is of weak order p = 1. If also the equations 

2. a T Al=±, 3. a J {b\ + bl) = \, 4. aJbi = \ 

are fulfilled, then the SRK method is of weak order p — 2. Finally, if additionally 

5. a J A 2 l = \, 6. a T (Al) 2 = ±, 7. a J A(b 2 + bl) = \, 

l 



8. a T (b l (Abi)+b 2 (Ab 2 )) = l 9. a T Ab l = z 

10. a T {{Al)(b 2 + b 2 )) = l 11. a T ((Al)foi) = i, 



12. a J (b\ + b 2 2 ) 2 = \, 13. a J (b\ + b\b 2 2 ) = \ , 14. a J b\ = \, 
15. (a T b 2 ) 2 = ± 

are fulfilled, then the SRK method is of weak order p = 3. 

Proof First, we note that E y/<p (u) =0 for all trees u € U" dd which have an odd num- 
ber of stochastic nodes of one color, see [6| or also [4, 16|. For those of these trees 
which have an order p(u) < p + 1, by construction of the method and due to the as- 
sumptions on Jk, k — I,..., 2m, it holds also E y/0 (u) =0. Thus, in the following we 
only have to consider trees with an even number of each kind of stochastic nodes, in 
particular only trees of integer order. Consequently, there are only two kinds of trees 
of order one to consider: 



u\= \> /, j = l,...,m, and u 2 




Theorems |2] and |3] O and ([33j yield 

V<p(u l )=hJ 2 , y l!> {u\)=Wj(h) 2 , y<j>(u 2 ) = ha T l and \y (!? {u 2 )=h. 

Thus, by the assumptions on Jj, E y/<p (u\ ) (h) =E (ui) (h) is fulfilled automatically, 
whereas E\]f$(u 2 )(h) = E\j/^,(u 2 )(h) yields order condition 1. 

If u £ Uj dd with u — [x\, . . . , T K ]f can be split into two trees u\ — [t ;i , . . . , T,-^ ]/, 
u 2 = [fj, ,Tj K ]/ with disjoint stochastic nodes, i. e. such that Kj , K 2 > 0, K\ + k 2 — 
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Table 3.1 Relevant second order trees and derivation of corresponding order conditions 



u 


v*(») 


V«>(») 
E VV(") 


ord. cond. 


(T) (T) 


CO XL) 
j 


h 2 J* 


Wj{hf 
3h 2 


by assumption 


C 

c 


7) 


Ira 1 At 
h 2 a T At 


tfsds 

h- 

2 


2. 




h 2 a T (b l J J + b 2 J m+] ) 2 
h 2 a J {b] + b 2 ) 


5liWj(.s) 2 ds 

lr 

2 


3. 


1 


; 


h 2 J j a T (b l Jj + b 2 J m+ j) 
h 2 a T bi 


Wj(h)jSWj(s)ds 
h- 

2 


4. 



fc, {j'i, . . . , i Kl ,j\ , . . . , j K2 } = {1, • • • , K - }, and the sets of colors of the stochastic nodes 
of u\ and ui are disjoint, then 

Ei^$(h)(/i) =Eip0(Mi)(/2)Ey#(K2)W =Ey< ! >(Ki)(/z) E Vp("2)v , = Ei/fy («)(/?), 

provided that the order conditions of orders lower than p (u) are fulfilled. Thus, in the 
following we only have to consider trees of second and third order which cannot be 
decomposed into two trees with disjoint stochastic nodes. The relevant second order 
trees together with the derivation of the corresponding order conditions are given in 
Table 13. ll the ones of order three in Tables [3~l2all3.2cl which completes the proof. 

Possible discrete choices for the random variables Jj c ,k= 1 , . . . , 2m, can be found in 
Table [331 



4 A concrete explicit third order SRK method 

Based on Theorem|4] we now calculate the coefficients of an explicit third order SRK 
method. The coefficients will be arranged in an extended Butcher array of the form 



c 


A 


h 


b 2 




a T 
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Table 3.2a Relevant third order trees and derivation of corresponding order conditions, part 1 



u 


v*(«) 




ord. cond. 




h 3 j) 
15/1 


I j/2 


by assumption 


® 








(a) 
© 


/i 3 a T A 2 l 
h 3 a T A 2 l 


Jo Jo' * 2 ds2 dsi 

A 3 
6 


5. 












h 3 a T (At) 2 
h 3 a T (At) 2 


h 3 
3 


6. 


^^^^ 


tfa T A{biJj + b2J m+j f 
h 3 a J A{b\+bl) 


;, 3 

6 


7. 



Whereas in the deterministic case we would only need three stages to construct an 
explicit third order method, here we need four stages to fulfill the 15 order conditions 
of Theorem[4] Therefore, we consider s = 4 in il.2\ , but require in addition that the 
method fulfills also the deterministic order four conditions. The remaining degrees 
of freedom are then eliminated by minimizing the vector lec of the order four coeffi- 
cients of the local error (13. Il l in the Euclidean norm assuming two dimensional noise 
(m = 2), i. e. by minimizing ||/ec||2 where 



lec = (P(u)-E[Y4>(u)(h) - V(p( u )( h )]) ue uf d .p(u)=4 ■ 



Using again the B-series analysis, a tedious calculation (one obtains 52 non automat- 
ically vanishing terms) and a subsequent attempt of numerical optimization yield the 
scheme AN3D1 presented in Table 14. ll AN3D1 needs two random variable and four 
drift evaluations per step, and thus two random variable and three drift evaluations 
less than Platen's third order method. 
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Table 3.2b Relevant third order trees and derivation of corresponding order conditions, part 2 



u 


Ev*(») 


<M") 
Ew(h) 


ord. 
cond. 


© @ 

I 


A 3 a T ( (6 1 Jj + b 2 J m+ j)(A(biJj+b 2 J m +j))) 
h i a J {b l {Ab x )+b 1 (Ab 2 )) 




8. 


® 

© © 


h 3 Jj(X T A {b\Jj + b 2 J m + j ) 
h 3 a T Ab! 


}? 

6 


9. 


©^^3^CD 


^^{(Arii^Jj + hJm+j) 2 ) 
^a T ({Al)(bj +bj)) 


JosWj(i)*ds 
;, 3 

3 


10. 




/iV,a T ((Al)(fe 1 7 ; + fo 2 7,„ + ,)) 
/i 3 a T ((Al)fo 1 ) 


3 


11. 




/i 3 a T ((fti/ ; +& 2 7 m+ ;) 2 (ii4+M n+ *) 2 ) 
\h i a T (?,b\+6b\bl + ?,b\) j = k 
j/i 3 a T (2> 2 +& 2 ) 2 


Jjwj(*) a wi(*) 2 «fe 

If j** 


12. 



5 Numerical example 

In the following we compare for three simple test equations the performance of the 
SRK scheme AN3D1 (with N(0,l)-distributed random variables) presented in the last 
section with some well known schemes, namely the third and the second order SRK 
schemes due to Platen 1 1 1|, denoted here by PL3 and PL2, respectively, DRI1 due to 
Debrabant and RoBler [9], and the extrapolated Euler-Maruyama scheme EXEM (cp. 
|20|) also attaining order two, which is given by 2E( {Z h l 2 (f )) ) - E( (Z h (f )) ), based 
on the Euler-Maruyama approximations Z h l 2 (t) and Z h (t) calculated with step sizes 
h/2 and h. In each case, the functional u(t) = E(f(X(t))) is approximated by a Monte 
Carlo simulation. The sample average mm,a(0 = jj Litli / (Y h (t, (Ok)), (Ok £ of 
M = 10 9 independent simulated realizations of the considered approximation Y h (t) 
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Table 3.2c Relevant third order trees and derivation of corresponding order conditions, part . 



u 


V*(") 




ord. 
cond. 


w 


h 3 J k a T ((b 1 Jj + b 2 J m+j ) 2 (b 1 J k + b 2 J m+k )) 
f/i 3 a T (3fcf+3&i&l) j = k 
\h 3 a r (b{+bibl) jjtk 


W k (h)J*Wj(s) 2 W k (s)ds 
jh 3 j = k 

U it* 


13. 




h 3 J j J k a T ((biJj+b 2 J m+ j)(b 1 J k + b 2 J m+k ) ) 
(h 3 a T (3b 2 +b 2 2 ) j = k 
\h 3 a T b 2 j^k 


Wj(h)W k (h)jSWj{s)W k (s) ds 


due 
to 3. 

equiv. 

to 14. 




/i 3 ;]a T (i)i/ ; +W m+J ) 
3h 3 a T b l 


Wj(h) 3 f»Wj(s)ds 

3ft 3 
2 


4. 




h 3 {a T (biJj + b 2 J m+J )) 2 
/i 3 ((a T fci) 2 + (a T & 2 ) 2 ) 


ft 3 

3 


due 
to 4. 

equiv. 

to 15. 



Table 3.3 Some discrete random variables corresponding up to the ith moment to N(0, 1) 



i\distribution 



1 P(J k = 0) = 1 

?, P{J k = l)=P{J k = -\)= l 1 

5 P(J k = V3) = />(/* = -V3) = i P(A = 0) = \ 

7 P(J k = V6)=P(J k = -V6) = ±, p(J k = t)=P(J k = -l) = ±, P(7, = 0) = i 



Table 4.1 Coefficients of AN3D1 






















1 


1 













b 2 


1/2 


3/8 


1/8 








hi 


1 


-0.4526683126055039 


-0.4842227708685013 


1.9368910834740051 











1/6 


-0.005430430675258792 


2/3 


0.1720970973419255 





with 



bi = (-0.01844540496323970,0.8017012756521233. 0.5092227024816198, 0.9758794209767762) T 
b 2 = (-0.1866426386543421, -0.8575745885712401, -0.4723392695015512,0.3060354860326548) T 
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Fig. 5.1 Computational effort per simulation path versus precision for SDE j5.lt 

is calculated in order to estimate the expectation and thus to determine the systematic 
error of the considered schemes. In the following, we denote by fx = umm (T) — u(T) 
the mean error at time T and by 6^ the empirical variance of the mean error. Further, 
we calculate the confidence interval with boundaries a and b to the level of 90% for 
the estimated error p. (see [ 1 1 1 for details). 

First, we compute the second moment of the solution of the linear SDE 



The solution value E(X 2 (T)) is now approximated with step sizes 2 , ... ,2~ 4 at time 
T = 2. The results for the applied schemes are presented in Table 15.11 Of course, 
these results have to be related to the computational effort of the schemes which we 
take in the following as sum of the number of evaluations of the drift function a 
as well as the number of random variables that have to be simulated. Then we can 
oppose the computational efforts to the errors of the analyzed schemes. The results 
are presented in Figure lBTTI Although being of different order, the two Platen schemes 
yield comparable results. This is due to the much higher computational costs of PL3. 
Both methods are better than the extrapolated Euler method, but perform worse than 
DRI1, which has optimized coefficients [9] and behaves therefore nearly like an order 
three method. Our new method AN3D1 performs best. 
As next example we consider the nonlinear SDE 




(5.1) 



which can be calculated analytically as 





(5.3) 
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Table 5.1 Mean errors, empirical variances and confidence intervals for SDE )5.U 





h 


A 




a 


b 




2 1 


-1.900E+02 


7.882E-07 


-1.900E+02 


-1.900E+02 


EXEM 


2° 


-1.499E+02 


7.032E-06 


-1.499E+02 


-1.499E+02 


2" 1 


-9.357E+01 


3.666E-05 


-9.357E+01 


-9.357E+01 




2- 2 


-4.435E+01 


8.881E-05 


-4.435E+01 


-4.434E+01 




2- 3 


-1.649E+01 


1.988E-04 


-1.650E+01 


-1.649E+01 




2" 4 


-5.170E+00 


2.770E-04 


-5.174E+00 


-5.166E+00 




2 1 


-1.840E+02 


6.783E-07 


-1.840E+02 


-1.840E+02 


PL2 


2° 


-1.294E+02 


6.348E-06 


-1.294E+02 


-1.294E+02 


2-' 


-6.412E+01 


2.790E-05 


-6.412E+01 


-6.412E+01 




2- 2 


-2.312E+01 


4.995E-05 


-2.312E+01 


-2.312E+01 




2- 3 


-6.880E+00 


5.861E-05 


-6.882E+00 


-6.878E+00 




2" 4 


-1.863E+00 


8.264E-05 


-1.865E+00 


-1.861E+0O 




2 1 


-8.377E+01 


2.936E-03 


-8.378E+01 


-8.375E+01 


PL3 


2° 


-2.705E+01 


1.614E-03 


-2.706E+01 


-2.704E+01 


2" 1 


-8.941E+00 


2.345E-04 


-8.944E+00 


-8.937E+00 




2- 2 


-1.951E+00 


6.624E-05 


-1.953E+00 


-1.949E+00 




2" 3 


-3.111E-01 


6.180E-05 


-3.130E-01 


-3.093E-01 




2 -4 


-4.307E-02 


4.718E-05 


-4.470E-02 


-4.144E-02 




2 1 


-1.316E+02 


3.486E-06 


-1.316E+02 


-1.316E+02 


DRIl 


2 U 


-5.438E+01 


2.049E-05 


-5.438E+01 


-5.437E+01 


2-' 


-1.308E+01 


4.872E-05 


-1.308E+01 


-1.308E+01 




2- 2 


-2.254E+00 


6.097E-05 


-2.256E+00 


-2.252E+00 




2- 3 


-3.314E-01 


6.225E-05 


-3.333E-01 


-3.295E-01 




2 -4 


-4.343E-02 


8.405E-05 


-4.560E-02 


-4.126E-02 




2 1 


-7.638E+01 


1.286E-05 


-7.638E+01 


-7.638E+01 


AN3D1 


2° 


-1.654E+01 


4.729E-05 


-1.654E+01 


-1.654E+01 


2" 1 


-1.946E+00 


6.804E-05 


-1.948E+00 


-1.944E+00 




2-2 


-1.651E-01 


4.993E-05 


-1.668E-01 


-1.635E-01 




2- 3 


-1.073E-02 


5.940E-05 


-1.255E-02 


-8.900E-03 




2-4 


-1.030E-04 


4.754E-05 


-1.738E-03 


1.532E-03 



Then E(e 2jr W) can be calculated as 

E(« K «)) = (V/5 + ™\ e ioi/50i _ 150_ (5 4) 

y " \ ioi y 101 

The solution value E(e 2x ' r ') is approximated with step sizes 2 1 , . . . ,2~ 4 at time T = 
2. The results for the applied schemes are presented in Table 15.21 and Figure [5721 and 
reflect a similar behaviour to the one from the linear example, except that PL3 suffers 
now from stability problems. 

As last example we consider the following linear system of SDEs with two di- 
mensional noise 

^-GK(4A)™* + tti)» —a* 



(5.5) 



where Epfj (f)) can be calculated as 



, 2 37 + 31148e- 5 '/ 4 -1185e-' 



E(Xf(t)) = ■ . (5.6) 

v 2 K " 30000 
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Table 5.2 Mean errors, empirical variances and confidence intervals for SDE 15.31 





h 


A 




a 


b 




2 1 


-7.925E+03 


2.818E-01 


-7.925E+03 


-7.925E+03 


EXEM 


2° 


-4.127E+02 


2.748E-03 


-4.127E+02 


-4.127E+02 


2" 1 


-4.777E+01 


8.555E-04 


-4.777E+01 


-4.776E+01 




2- 2 


-7.296E+00 


6.418E-04 


-7.302E+00 


-7.290E+00 




2- 3 


-1.369E+00 


5.752E-04 


-1.375E+00 


-1.363E+00 




2" 4 


-2.992E-01 


5.557E-04 


-3.048E-01 


-2.936E-01 




2 1 


6.573E+02 


2.146E-03 


6.573E+02 


6.573E+02 


PL2 


2° 


1.010E+02 


1.921E-04 


1.010E+02 


1.011E+02 


2-' 


1.678E+01 


1.192E-04 


1.678E+01 


1.678E+01 




2- 2 


2.676E+00 


1.064E-04 


2.674E+00 


2.679E+00 




2- 3 


4.665E-01 


1.064E-04 


4.641E-01 


4.690E-01 




2" 4 


9.702E-02 


1.034E-04 


9.461E-02 


9.943E-02 




2 1 


8.202E+10 


1.589E+23 


-1.249E+10 


1.765E+11 


PL3 


2° 


Inf 


NaN 


NaN 


NaN 


2" 1 


NaN 


NaN 


NaN 


NaN 




2- 2 


3.810E-01 


6.677E-05 


3.791E-01 


3.829E-01 




2- 3 


-9.317E-02 


8.318E-05 


-9.533E-02 


-9.101E-02 




2" 4 


-1.930E-02 


9.220E-05 


-2.158E-02 


-1.703E-02 




2 1 


1.360E+03 


9.350E-02 


1.360E+03 


1.360E+03 


DRIl 


2 U 


3.948E+01 


6.535E-05 


3.947E+01 


3.948E+01 


2-' 


1.525E+00 


8.976E-05 


1.522E+00 


1.527E+00 




2-2 


-1.412E-01 


1.012E-04 


-1.436E-01 


-1.388E-01 




2- 3 


-3.861E-02 


1.054E-04 


-4.105E-02 


-3.618E-02 




2 -4 


-3.432E-03 


1.032E-04 


-5.841E-03 


-1.022E-03 




2 1 


3.649E+01 


8.416E-05 


3.648E+01 


3.649E+01 


AN3D1 


2° 


1.871E+00 


8.809E-05 


1.869E+00 


1.873E+00 


2" 1 


-4.186E-01 


7.102E-05 


-4.206E-01 


-4.166E-01 




2-2 


-6.042E-02 


6.130E-05 


-6.228E-02 


-5.857E-02 




2-3 


-5.103E-O3 


8.299E-05 


-7.263E-03 


-2.943E-03 




2 -4 


3.022E-06 


9.237E-05 


-2.276E-03 


2.282E-03 











-«-EXEM 




- - 9-PL2 




...... PL3 




- -«- DRI1 




+ AN3D1 





' 6 0.5 1 1.5 2 2.5 3 

log ( function and random number evaluations ) 



Fig. 5.2 Computational effort per simulation path versus precision for SDE j5.3t 
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Table 5.3 Mean errors, empirical variances and confidence intervals for SDE )5.5t 





h 


A 


< 


a 


b 




2 1 


-3.122E-01 


1.279E-10 


-3.122E-01 


-3.121E-01 


EXEM 


2° 


-7.717E-03 


1.808E-11 


-7.718E-03 


-7.716E-03 


2-' 


-1.032E-03 


2.299E-11 


-1.033E-03 


-1.030E-03 




2- 2 


-1.848E-04 


2.527E-11 


-1.860E-04 


-1.836E-04 




2 -3 


-3.759E-05 


2.876E-11 


-3.887E-05 


-3.632E-05 




2 1 


3.491E-01 


9.822E-12 


3.491E-01 


3.491E-01 


PL2 


2° 


2.984E-02 


I.040E-11 


2.984E-02 


2.984E-02 


2" 1 


4.796E-03 


8.046E-12 


4.796E-03 


4.797E-03 




2- 2 


1.006E-03 


7.765E-12 


1.005E-03 


1.006E-03 




2 


2.325E-04 


6.647E-12 


2.319E-04 


2.331E-04 




2 1 


-4.468E-02 


3.591E-13 


-4.468E-02 


-4.468E-02 


DRIl 


2° 


-4.712E-03 


5.216E-12 


-4.713E-03 


-4.712E-03 


2- 1 


-4.603E-04 


7.853E-12 


-4.610E-04 


-4.597E-04 




2- 2 


-5.199E-05 


7.022E-12 


-5.262E-05 


-5.137E-05 




2- 3 


-6.531E-06 


6.369E-12 


-7.130E-06 


-5.933E-06 




2' 


2.526E-02 


9.516E-12 


2.526E-02 


2.526E-02 


AN3D1 


2° 


7.390E-04 


7.190E-12 


7.383E-04 


7.396E-04 


2- 1 


3.150E-05 


5.109E-12 


3.096E-05 


3.204E-05 




2- 2 


1.459E-06 


6.996E-12 


8.314E-07 


2.086E-06 




2- 3 


4.859E-08 


6.152E-12 


-5.395E-07 


6.367E-07 



^ -3 





















***** 




-«-EXEM 








-3-PL2 








-«- DRI1 








+ AN3D1 









8.5 1 1^5 2 2.5 

log ( function and random number evaluations ) 



Fig. 5.3 Computational effort per simulation path versus precision for SDE j5.5t 



The solution value ^{X^T)) is approximated with step sizes 2',...,2~ 3 at time 
T = 2. The results are presented in Table 1531 and Figure [33] (note that PL3 is not 
applicable here). Again, AN3D1 performs best. 



6 Conclusion 

We have presented a general class of SRK methods for the weak approximation of 
SDEs with additive noise, together with the corresponding order conditions up to 
order three. A concrete explicit third order method has been derived, for which a nu- 
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merical comparison with some well known other methods regarding its performance 
yielded very promising results. In contrast to the method of Platen, it needs only two 
random variables and four drift evaluations per step and is also applicable to SDEs 
driven by a multidimensional Wiener process. Future research may be done by con- 
structing implicit methods with good stability properties, i. e. which are suitable for 
stiff problems, and by developing methods for more general noise. 

Acknowledgement 

The author is grateful to Birgit Debrabant and an anonymous referee for their helpful 
hints which improved the presentation of the material. 

References 

1. Burrage, K., Burrage, P.M.: High strong order explicit Runge-Kutta methods for stochastic ordinary 
differential equations. Appl. Numer. Math. 22(1-3), 81-101 (1996). Special issue celebrating the 
centenary of Runge-Kutta methods 

2. Burrage, K., Burrage, P.M.: Order conditions of stochastic Runge-Kutta methods by B-series. SIAM 
J. Numer. Anal. 38(5), 1626-1646 (electronic) (2000) 

3. Burrage, K., Burrage, P.M., Tian, T.: Numerical methods for strong solutions of stochastic differential 
equations: an overview. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 460(2041), 373^02 (2004). 
Stochastic analysis with applications to mathematical finance 

4. Burrage, P.M.: Runge-Kutta methods for stochastic differential equations. Ph.D. thesis, The Univer- 
sity of Queensland, Brisbane (1999) 

5. Debrabant, K., Kvsern0, A.: B-series analysis of stochastic Runge-Kutta methods that use an iterative 
scheme to compute their internal stage values. SIAM J. Numer. Anal. 47(1), 181-203 (2008/09) 

6. Debrabant, K, Kvsern0, A.: Stochastic Taylor expansions: Weight functions of B-series expressed as 
multiple integrals. Stoch. Anal. Appl. 28(2), 293 - 302 (2010). DOI 10.1080/07362990903546504 

7. Debrabant, K., RoBler, A.: Continuous weak approximation for stochastic differential equations. J. 
Comput. Appl. Math. 214(1), 259-273 (2008) 

8. Debrabant, K., RoBler, A.: Diagonally drift-implicit Runge-Kutta methods of weak order one and two 
for Ito SDEs and stability analysis. Appl. Numer. Math. 59(3-4), 595-607 (2009) 

9. Debrabant, K., RoBler, A.: Families of efficient second order Runge-Kutta methods for the weak 
approximation of Ito stochastic differential equations. Appl. Numer. Math. 59(3-4), 582-594 (2009) 

10. Karatzas, I., Shreve, S.E.: Brownian motion and stochastic calculus, Graduate Texts in Mathematics, 
vol. 113, second edn. Springer- Verlag, New York (1991) 

11. Kloeden, P.E., Platen, E.: Numerical solution of stochastic differential equations, Applications of 
Mathematics, vol. 21, 2 edn. Springer- Verlag, Berlin (1999) 

12. Komori, Y.: Multi-colored rooted tree analysis of the weak order conditions of a stochastic Runge- 
Kutta family. Appl. Numer. Math. 57(2), 147-165 (2007) 

13. Komori, Y: Weak second-order stochastic Runge-Kutta methods for non-commutative stochastic dif- 
ferential equations. J. Comput. Appl. Math. 206(1), 158-173 (2007) 

14. Mackevicius, V., Navikas, J.: Second order weak Runge-Kutta type methods of Ito equations. Math. 
Comput. Simulation 57(1-2), 29-34 (2001) 

15. Milstein, G.N.: Numerical integration of stochastic differential equations, Mathematics and its Ap- 
plications, vol. 313. Kluwer Academic Publishers Group, Dordrecht (1995). Translated and revised 
from the 1988 Russian original 

16. RoBler, A.: Stochastic Taylor expansions for the expectation of functionals of diffusion processes. 
Stoch. Anal. Appl. 22(6), 1553-1576 (2004) 

17. RoBler, A.: Rooted tree analysis for order conditions of stochastic Runge-Kutta methods for the weak 
approximation of stochastic differential equations. Stoch. Anal. Appl. 24(1), 97-134 (2006) 

18. RoBler, A.: Second order Runge-Kutta methods for Stratonovich stochastic differential equations. 
BIT 47(3), 657-680 (2007) 



16 



Kristian Debrabant 



19. RoBler, A.: Second order Runge-Kutta methods for Ito stochastic differential equations. SIAM J. 
Numer. Anal. 47(3), 1713-1738 (electronic) (2009) 

20. Talay, D., Tubaro, L.: Expansion of the global error for numerical schemes solving stochastic differ- 
ential equations. Stoch. Anal. Appl. 8(4), 94-120 (1990) 

21. Tocino, A., Vigo-Aguiar, J.: Weak second order conditions for stochastic Runge-Kutta methods. 
SIAM J. Sci. Comput. 24(2), 507-523 (electronic) (2002) 



